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Abstract 

The F 4- H 2 — ► HF + H potential energy hypersurface has been studied in the 
saddle point and entrance channel regions. Using a large [5s 5p 3 d 2/ 1 g/4s 3 p 2d] 
atomic natural orbital basis set, we obtain a classical barrier height of 1.86 kcal/mole 
at the CASSCF /multireference Cl level (MRCI) after correcting for basis set su- 
perposition error and including a Davidson correction (-f Q) for higher excitations. 
Based upon an analysis of the computed results, the true classical barrier is esti- 
mated to be about 1.4 kcal/mole. We also compute the location of the bottleneck 
on the lowest vibrationally adiabatic potential curve, and determine the transla- 
tional energy threshold from a one-dimensional tunneling calculation. Using the 
difference between the calculated and experimental threshold to adjust the classi- 
cal barrier height on the computed surface yields a classical barrier in the range 
of 1.0-1. 5 kcal/mole. Combining the results of our direct estimates of the classical 
barrier height with the empirical values obtained from our approximate calcula- 
tions of the dynamical threshold, we predict that the true classical barrier height 
is 1.4± 0.4 kcal/mole. Arguments are presented in favor of including the relatively 
large (ssl kcal/mole) +Q correction obtained when nine electrons are correlated at 
the CASSCF /MRCI level. 
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I. Introduction 


The F + H 2 -* HF + H reaction is a prototype in the field of gas phase re- 
action dynamics, since it is accessible to both detailed theoretical and experimen- 
tal study. Early electronic structure studies of this reaction by Bender et al. [l], 
using first-order Cl wave functions in a DZP basis set, gave a barrier height of 
1.7 kcal/mole, but recent more extensive ab initio calculations by Frisch tt al. [2] 
gave a barrier height greater than 3.2 kcal/mole. In a recent review article, Schae- 
fer [3] reiterated his position that the energy barrier for this reaction on the low- 
est Born-Oppenheimer potential energy surface (PES) must be greater than the 
2.35 kcal/mole value deduced in an earlier study by Ungemach et al. [4]. Further, us- 
ing a harmonic vibrational analysis, it was argued that this “classical” barrier height 
is consistent with the observed experimental translational energy threshold [5] for 
the reaction (about 0.8 kcal/mole) when zero-point corrections are included. On 
the other hand, Steckler et al. [6] carried out canonical variational transition state 
theory (VTST) calculations, including the effects of vibrational anharmonicity and 
multidimensional tunneling corrections, using a global empirical potential energy 
surface with a 2.7 kcal/mole barrier and other saddle point properties based on 
the Frisch et al. [2] surface. They concluded that the barrier from that work must 
be significantly too high, since the VTST calculations produced a thermal rate 
constant at 190K that was about 30 times smaller than experiment. They also 
argued that the zero-point vibrational energy correction is much smaller, because 
the VTST bottleneck for the reaction is shifted much further out in the entrance 
channel from the classical barrier, and showed that it is important to include anhar- 
monic effects in the threshold determination. Recently, Truhlar and coworkers [7-9] 
have reexamined the saddle point region on the F+H 2 PES using a multireference 
Cl (MRCI) description. They attempt to reach both the one-particle and n-particle 
limits using the scaled external correlation (SEC) method [10]. For collinear geome- 
tries, their SEC energy barrier is 1.6 kcal/mole, much lower than previous MRCI 
calculations as well as recent studies using the quantum Monte Carlo method [ll]. 
Also, they find that the MRCI-f Q (where -f Q denotes the multireference Davidson 
correction for higher excitations) and the SEC bending potentials exhibited addi- 
tional minima for nonlinear geometries. Their extremely fiat bending potential is 
significantly different from those implicit in empirical LEPS surfaces such as the 
commonly used Muckerman surface No. 5 [12], and the first-order Cl calculations 
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of Bender et ai [13]. 

In this work we consider several aspects of the F + H 2 PES. First, we present 
a study of the classical barrier height as a function of improvements to both the 
one-particle and n-particle treatments. Recent advances [14] in strategies for basis 
set contraction, which allow the use of very large primitive basis sets, enable us 
to approach the one-particle basis set limit. We also are able to use a higher 
level of correlation for calibration than in previous treatments. Second, using the 
externally contracted Cl method (CCI) [15], bending potentials were computed in 
the collinear saddle point region. Our best computed bending potential favors a 
non-collinear approach, but by a lesser extent than does the SEC bending potential 
of Schwenke et al. [8]. Third, the calculated CCI surface was used to locate the 
bottleneck on the vibrationally adiabatic potential curve, and the reaction threshold 
was deduced from a one-dimensional tunneling calculation. We estimate the true 
classical barrier height by adjusting the CCI barrier height for the difference in 
the calculated and experimental thresholds. In the next section we discuss the 
details of the calculations such as basis sets, choice of active spaces in the orbital 
determinations and correlation treatments. In Section III we discuss our study of 
the one-particle and n-particle requirements for a direct ab initio computation of 
the classical barrier height. In Section IV we discuss our attempts to characterize 
an accurate vibrationally adiabatic potential curve for both F + H 2 and F + D 2 . 
Our conclusions are presented in Section V. 

II. Computational Details. 

The primitive basis for fluorine is the (13s 8p) set of van Duijneveldt {16], aug- 
mented with a (6 d 4/ 2 g) polarization set. The polarization functions are taken as 
even-tempered sequences with an internal ratio of 2.5. The geometric mean of the d 
exponents is 1.62, while the / and g sets are based on multiplying this mean d expo- 
nent by 1.2 and 1.44, respectively. This primitive set is contracted to [5s 4p 3d 2/ Iff] 
using a general contraction scheme based on atomic natural orbitals (AN Os) [14]. 
This basis is further augmented by an additional even-tempered diffuse 2 p function 
to better describe F~ character (o=0.059326). The primitive basis for hydrogen 
is the (8s) set of van Duijneveldt [16], augmented with a (6p 4d) polarization set. 
The p and d sets are even-tempered sequences with an internal ratio of 2.5 and 
geometric means of 1.0. The hydrogen / basis set consists of three even-tempered 
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primitive / functions with a geometric mean exponent of 1.0. This primitive set 
is contracted to [4s 3 p 2d If] using ANOs from a calculation on H 2 , as described 
by Almlof and Taylor [14]. The sets described here are the largest contracted sets 
used in the present work: smaller sets were obtained by deleting ANOs as required. 
The 3s, 4p, 5s and 5 d combinations of the 3d, 4/ and 5 g functions, respectively, are 
deleted in all calculations. 

For all collinear geometries, the SCF configuration of the F+H 2 2 E + surface 
connecting reactants and products is: 

la 7 2a 7 3a 7 4a 1 In 4 , 

where la and 2 a are nominally the fluorine Is and 2s orbitals. In the reactants 
channel, the 3a orbital is the H 2 bonding orbital (la c ), and the 4a and lx orbitals 
are predominantly the fluorine 2p orbitals. This reference configuration was used 
for all of our single-reference based correlation treatments, which included singles 
plus doubles Cl (SDCI) and the coupled pair functional (CPF) method [17]. 

CASSCF wave functions were used as the starting point for our multireference 
Cl (MRCI) treatments. The minimum active space that gives a reasonable first- 
order description of this system comprises the 3o-5a orbitals. For the reactants 
channel, the 5a orbital corresponds to the H 2 la u orbital. In C 2 V symmetry this is 
denoted as a (300) active space, where the integers denote the number of active ai, 
bi and b 2 orbitals. When the lx orbital is included in the active space to give a (311) 
CASSCF, the number of configurations is increased only slightly. The barrier is not 
significantly changed when the MRCI calculation is based on the larger (311) space. 
It is well known that 2p — ♦ 2 p' correlation is important in accurately computing the 
electron affinity (EA) of F and O [18], as well as in describing HF and OH owing 
to the large fraction of F - and O - character in the wave function. In the course 
of the reaction X+H 2 — >HX +H, the degree of X - character in the wave function 
increases. It is, therefore, not surprising that Liu [19] has found that inclusion 
of 2p — » 2p' excitations in the CASSCF and MRCI reference space significantly 
(0.8 kcal/mole) reduced the barrier height for the F+H 2 reaction. This has also 
been found for 0( 3 P)+H2 by Walch [20]. Thus to accurately compute the F+H 2 
barrier, the 2p' shell must be included in the CASSCF resulting in a (422) active 
space. However, the fourth a } orbital is only weakly occupied and essentially of the 
same importance as the 2 a g orbital in H 2 . Thus either a (322) or (522) CASSCF 
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active space must be used to eliminate the possibility of multiple solutions. The 
(522) active space is sufficiently flexible to correlate both the F 2 p orbital and to 
provide radial correlation for H 2 . Alternatively, as discussed below, a (322) active 
space, which takes into account most of the effect of 2p— »2p' excitations, can be 
used at less expense than the (522) active space. 

The calculation of a smooth bending potential about the collinear saddle point 
required an additional expansion of the active space. Two problems occurred when 
the (322) active space (denoted (52) in C e symmetry) was used for non-collinear 
geometries. First, this active space includes the 2p' orbitals, which correlate the 
doubly occupied fluorine 2p orbitals, but not the singly occupied one. For collinear 
geometries this is a reasonable choice based on natural orbital occupations, but for 
large bending angles this wave function changes discontinuously. For angles 6 > 
75°, (0=180° — Z(HHF)), the 2p a correlation is more important than the H 2 left- 
right correlation, which leads to a change in character of the most weakly occupied 
CASSCF orbital. Thus, the calculation using the col linear-derived active space 
is not adequate for larger 0 values needed to establish the shape of the potential 
curve. Second, there is a competition between 2s and 2 p correlation, which leads 
to a mixing of the fluorine 2s and the in-plane doubly occupied 2p orbital. This 
causes symmetry breaking of the orbitals for collinear geometries when treated 
in C, symmetry, resulting in a discontinuity between bent and linear structures. 
Adding the other component of 2 p' leads to a (62) active space, while also adding 
the fluorine 2s and 2s' orbitals leads to an (82) active space. As shown in Fig. 
1, smooth bending potentials are obtained only when the correlation treatment is 
based on the (82) active space — see Section IV. 

Multireference treatments using both the MRCI method and the externally con- 
tracted Cl (CCI) method of Siegbahn [15] have been performed from several active 
spaces. For the (322) and larger reference spaces, the complete second-order treat- 
ment is prohibitively large, thus necessitating reference selection. Therefore, we 
include an occupation in the reference space only if the absolute value of the coeffi- 
cient of one of its component spin couplings in the CASSCF wave function is above 
a designated threshold. Thresholds of 0.05, 0.025 and 0.01 are used in this work. 
In each case the multireference analog of Davidson’s correction [21 j was added to 
account for the contribution from unlinked quadruple excitations. The correction 
in the MRCI case is AE (1 — X^C^) where the Cr are the coefficients of the ref- 
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erence configurations in the MRCI wave function and AE is the energy difference 
between the reference and MRCI wave functions. In the CCI case this expression is 
divided by The results where Davidson’s correction has been applied are 

denoted by MRCI + Q or CCI + Q. 

Since the CCI method is less expensive than the MRCI method, it allowed us 
to use larger configuration state function (CSF) expansions, to consider a larger 
number of geometries on the F + H 2 surface, and to investigate the effects of further 
basis set extensions. At the CCI(CCI-fQ) level using a 0.025 selection threshold 
and the [5s 5p 3d 2/ 1 g/As 3 p 2d] ANO basis set, the collinear barrier height is 
2.79(2.02) kcal/mole, which is about 0.2— 0.4 kcal/mole larger than at the MRCI 
level. Removing the g function leads to a barrier height of 2.89(2.15) kcal/mole. 
This 0.1 kcal/mole basis set expansion error was deemed small enough that the 
[5s 5p 3d 2/ /4s 3 p 2d] basis set was used to generate a 20-point grid in rjjjr and Tgg, 
which spanned the saddle point and dynamical bottleneck for the entrance channel 
region of the PES. To verify that the restrictions inherent in the CCI approach do 
not significantly alter the shape of the surface, MRCI calculations were carried out 
for nine collinear geometries centered around the saddle point and at the bottleneck 
geometry. The MRCI+Q and CCI+Q quadratic force constants at the saddle point 
differ by less than 10%. The difference in energy between the saddle point and 
bottleneck geometry is 0.19 kcal/mol (MRCI+Q) versus 0.28 kcal/mole (CCI+Q). 
The effect of these uncertainties in the CCI+Q potential on the computed threshold 
is discussed in Section IV. 

The CCI potential calculations were carried out on the CRAY X-MP/48 with the 
MOLECULE [22]-SWEDEN [23] system of programs. The very large MRCI calibra- 
tion calculations were done primarily on the NAS CRAY 2. The CPF calculations 
were done on the Cyber 205 using the Karlsruhe codes [24]. 

III. Studies of the classical barrier height 

In Table I we compare, at the 7-electron level (His and F2p) of correlation, vari- 
ous MRCI (and MRCI+Q) treatments with the full configuration-interaction(FCI) 
treatment using the [4s 3p ld/2s lp] ANO basis set: this is an expansion of the 
approximate treatments of the correlation problem reported in Ref. 25. The entries 
are denoted both by the active space and by the reference threshold selection crite- 
rion. Most noteworthy is that the barrier heights are in much better agreement with 
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the FCI after a correction is made for quadruple excitations. Also, the magnitude 
of the +Q correction decreases as the reference threshold is lowered. Further, the 
barrier height varies much less with selection threshold if a -)-Q correction is in- 
cluded. For the larger (522) active space the -fQ correction does slightly overshoot 
the FCI result. Similar conclusions can be drawn from part B of Table I, where 
we have optimized the geometry at each level of treatment using a biquadratic fit 
to a grid of nine points. Both the geometry and barrier height are generally in 
better agreement with the FCI after the +Q correction is added. Although the -l-Q 
correction may be too large in the limit of a SOCI calculation (i.e. no reference 
selection), we have observed that the MRCI+Q results are often superior to the 
MRCI if the reference configurations are selected. We also believe that the +Q 
correction is even less likely to be too large when 9 electrons are correlated, but 
unfortunately we are unable to carry out the necessary FCI calibration calculations 
for a realistic one-particle basis. However, the observations made in our study of 
the electron affinity [26] of oxygen are probably also valid here. It appears that 2s 
correlation increases the electron affinity of oxygen (by as much as 2.5 kcal/mole), 
but that a very high level of correlation treatment is required to fully account for 
this effect. Our work suggests that to obtain a quantitative barrier height for F+H 2 
it is not only necessary to explicitly include 2 p — > 2 j/ correlation, but also to fully 
account for differential effects from 2s correlation. 

Based upon our calibration calculations, the threshold for reference selection in 
the ANO basis set calculations is chosen as 0.0 (no selection) for the MRCI(300) 
reference space and 0.025 in all other cases. For example, in the MRCI(322) treat- 
ment, the final reference list is obtained by merging the important occupations for 
F+H 2 , HF-fH and the MRCI(300)-|-Q saddle point barrier geometry. This leads 
to 12 reference occupations (29 CSFs), and an MRCI expansion of 1 343 112 CSFs 
when 9 electrons are correlated. 

In Table II we have summarized our results for the reactants H 2 and HF as well 
as our calculated basis set superposition errors (SE). Our results for both diatomics 
are in very good agreement with experiment [27]. For HF, the errors are only 0.002- 
0.003 ao for r e and 0.05-0.06 eV for D e , for calculations using the (222) active space. 
The basis set SE have been computed using the counterpoise method [28] including 
the full space of the ghost basis. The basis set superposition errors for H 2 in the 
presence of the F basis are only about 0.01 kcal/mole. The basis set superposition 
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errors for F in the presence of the H 2 basis are also quite small, especially if the 2s 
electrons are not correlated. The superposition errors for F~ in the presence of the 
H 2 basis are much larger than those for F. However, Mulliken population analysis 
shows that in the saddle point region, the wave function possesses only about 5% 
F~ character. Hence it is appropriate to use the SE for F in Table II to correct 
our barrier heights; approximately 0.1 and 0.2 kcal/mole should be added to the 
7-electron and 9-electron barrier heights, respectively. The values in the tables were 
not corrected, since the remaining basis set incompleteness effects are probably of 
the same magnitude and of opposite sign. 

In Table III we summarize our theoretical studies of the classical saddle-point 
geometry and barrier for the F + H 2 reaction. Our MRCI(300) treatment in the 
[5s 5p 3d 2/ \gjAs 3 p 2d] ANO basis yields a barrier height and saddle-point ge- 
ometry in good agreement with the large Slater basis set calculations of Frisch et 
a!. [2]. This is especially true after the Slater calculations are corrected for their 
larger basis set superposition error (see Table II). Adding a +Q correction to the 
MRCI(300) result reduces the barrier by about 1.2 kcal/mole and shifts the saddle 
point by over 0.1 ao further out in the entrance channel. Another substantial re- 
duction of the barrier height occurs when the MRCI is based upon the (322) active 
space, which includes 2 p — ► 2 p' excitations; the effect is very similar in magnitude 
to that found by Liu [19]. The +Q for this active space has a smaller, but still 
very large, effect on the barrier height. If 2s correlation is not included, the bar- 
rier height is slightly higher and the magnitude of the differential -|-Q correction 
is smaller by a factor of about one half. The exothermicity of the reaction is also 
increased substantially, but interestingly the saddle-point geometry does not change 
significantly. Of the two single-reference based entries in Table III, the CPF method 
is in better agreement with the MRCI(322)+Q than SDCI+Q for both the saddle 
point geometry and barrier, although it is inferior to the MRCI-based methods for 
the exothermicity. It should also be noted that the addition of the +Q correction at 
the MRCI(322) level reduces the size consistency error, making the supermolecule 
exothermicity nearly equal to the value computed using the fragment dissociation 
energies. 

In the remainder of Table III we present some calibration calculations using the 
CCI method. We find that adding an f ANO to the hydrogen basis decreases the 
barrier by 0.06(0.07) kcal/mole at the CCI(CCI+Q) level. Furthermore, removing 
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the g function from fluorine increases the barrier by 0.11(0.13) kcal/mole. Appar- 
ently, then, basis set improvements decrease the barrier height. Considering that 
the SCF and CASSCF barriers are quite large, the barrier height tends to decrease 
monotonically when improvements are made in the level of correlation treatment 
that recover successively larger percentages of the differential correlation energy 
between the reactants and the saddle point. In the following discussion we attempt 
to assign upper bounds of varying rigor to our ab initio calculations. A very conser- 
vative upper limit of 2.52 kcal/mole is given by the 7-electron MRCI(322)— SE+Q 
calculation, since this MRCI is calibrated against FCI calculations, and we know 
that 2s correlation reduces the barrier. Since the 7-electron +Q correction is not 
likely to be an overestimate, based on the FCI calibration calculations, the value of 
2.26 kcal/mole obtained by the 9-electron MRCI(322)— SE +Q (7-electron) is also 
a conservative upper bound. Since the -(-Q correction is less likely to overshoot 
as more electrons are correlated, we expect that our 9-electron MRCI(322)— SE+Q 
result of 1.86 kcal/mole is a realistic estimate. However, our best estimate for the 
barrier includes a correction for remaining basis set incompleteness, as discussed 
above, and a larger +Q correction, which we feel is justified from studies of the 
electron affinities of oxygen and fluorine. Thus instead of correcting for basis set 
superposition error we assume that the basis set limit is actually 0.1 kcal/mole less 
than the directly computed value. A justifiable lower bound of 1.35 kcal/mole for 
the barrier height is obtained by using 120% of the +Q correction. Although this 
latter value is below the directly computed values, the extrapolations seem well jus- 
tified based on analogous studies of the electron affinity of oxygen, which indicate 
that the full effect of higher excitations and 2s correlation are difficult to incorporate 
into the MRCI calculations. These lower values are consistent with extrapolations 
of Truhlar and coworkers [7-9] using the SEC method. Further, they are consistent 
with estimates obtained using VTST and our ab initio bending potentials as dis- 
cussed in the next section. Based only upon our ab initio calculations, the classical 
barrier height is, therefore, estimated to be between 1.35 and 1.86 kcal/mole. 

IV. Studies of the vibrationally adiabatic potential curve 

It is well known that, for gas phase reactions, the “classical” barrier height on the 
PES is only approximately related to the experimental activation energy or thresh- 
old energy. Indeed, for F+H 2 previous studies [2, 3, 7-9] have shown that zero-point 
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vibrational energy effects may be as large as 50% of the computed barrier height. 
Therefore, to better characterize the saddle point and entrance channel regions of 
the potential energy surface, a grid of 20 collinear geometries was obtained at the 
CCI(322) level. These points encompassed the dynamical bottleneck region. A 
bending potential was then computed with thf and thh held fixed at approxi- 
mately their collinear saddle point values (r/*jr =2.90 ao and r jjjj =1.45 a 0 ). An 
explicit tabulation of this CCI potential energy grid is available from one of the au- 
thors (SPW). The bending potentials computed with the (52), (62) and (82) active 
spaces (denoted in C, symmetry) are illustrated in Fig. 1. Each of these CCI+Q 
bending potentials was generated using a reference space that included all CSFs 
with coefficients greater than 0.025 in the CASSCF wave functions (merged list for 
representative bending angles). As the F 4 -H 2 geometry is not included in the selec- 
tion of references, these calculations are not appropriate for computing the actual 
barrier height, and hence are used only to compute the bending potential. Only the 
(82) curve shows a smooth variation with angle. The effects of symmetry breaking 
are illustrated by the (62) bending curve. The connected points for the (62) bend- 
ing curve were all performed using C s symmetry and this potential is similar to 
the (82) potential, but not as smooth. The single point in Fig. 1 is for a collinear 
calculation in C 2 V symmetry. If 0 = 0° were to pass through this point, instead 
of the broken symmetry point, then a curve very similar to the (52) curve would 
be obtained. Interestingly, the orbitals from the (52) CASSCF calculation do not 
break symmetry for collinear geometries computed in C, symmetry. However, the 
non-collinear geometries do exhibit mixing of the fluorine 2s and in-plane doubly 
occupied 2p orbitals, and the maximum at 15° in the (52) potential is believed to be 
nonphysical. It is interesting to note that the upper curves in Fig. 1 resemble the 
MRCI+Q and SEC bending potentials of Schwenke ct al. [8] in that they exhibit 
minima at 0=0° as well as at larger angles. 

The (82) bending potential, shown in Fig. 1, is very flat up to 6 «60°, but rises 
fairly sharply for larger 8 values. For this choice of active space, the CCI and CCI 
+ Q bending potentials are very similar. The angle variation is only carried to 
8 = 105°, because for that angle the structure has approximately C 2 V symmetry 
with the F nearly equidistant from the two H atoms. It is interesting to note that the 
energy at that geometry is only 2.37 kcal/mole higher than at the collinear saddle 
point in the CCI calculations (1.81 kcal/mole higher in CCI-f-Q). This means that 
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the potential energy surface is rather isotropic with respect to the F atom approach 
to the H 2 molecule. The actual minimum in the CCI bending potential occurs for 
6 «60°, but the bent structure is only 0.04 kcal/mole below the collinear saddle 
point structure (0.12 kcal/mole lower at the CCI+Q level). This bending potential 
is qualitatively similar to the one obtained by the SEC method [8], except that the 
SEC method predicts a stronger energetic preference for non-collinear geometries 
(by about 0.3 kcal/mole) them does the present calculation and has an additional 
energy minimum at 180°. In fact, the SEC bending potential more closely resembles 
the (52) CCI+Q bending potential drawn in Fig. 1. The Bender tt al. [13] bending 
potential also predicts the minimum energy approach to be slightly noncollinear, 
but their potential rises somewhat more steeply for larger 6 values. 

Canonical variational transition state theory (VTST) [29] has been shown to 
provide a good estimate of rate constants for simple gas phase reactions. In this 
approach, the transition state is located at the free energy “bottleneck" between 
reactants and products, and the rate constant is computed from transition state 
theory with corrections for tunneling and recrossing effects. To locate the variational 
transition state, the triatomic vibrational frequencies along the minimum energy 
path for F + H 2 collisions must first be determined. 

The bending contribution was determined from the one-dimensional bending po- 
tential by numerical integration of the vibrational Schrodinger equation with fixed 
r hh and r jj F . This computation used a standard diatomic vibrational energy code 
with the bending angle 6 scaled by bf t H H • The bending potential energy curve 

was represented by a bicubic spline fit. The effective reduced mass was taken to be 

_ rogmjrrgffrj/jr 

Me// rn H r 2 HH + m F r 2 HF + m F (r HU + riiF ) 2 ' 

which comes from the scaled Wilson G-matrix element [30] for the bending coor- 
dinate of a linear triatomic molecule. To validate this approach, we performed a 
test calculation for a quadratic-quartic potential with known energy levels. Our 
approach obtained excellent agreement for the three lowest energy levels, and only 
small inaccuracies for the higher levels. The resulting energy levels for the F-H-H 
and F-D-D systems with the CCI+Q bending potential shown in Fig. 2 demonstrate 
that the energy levels are rather anharmonic, and that the zeroth level lies slightly 
above the energy maximum found in the bending potential at collinear geometry. 
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Therefore the reference reaction path for this reaction can be taken to be collinear 
even though the bending potential shows a slight preference for bent geometries. 
As the contribution to the zero-point energy from the bending mode in F + H 2 and 
F + D 2 at the collinear saddle point is fairly small, and goes to zero for the sep- 
arated reactants, it is a reasonable assumption to consider it to be constant along 
the minimum energy path (MEP) in the vicinity of the dynamical bottleneck. 

The next step in locating the variational transition state is to determine the 
symmetric and asymmetric stretching frequencies along the MEP. To do this we 
fitted the 4x5 grids of collinear CCI and CCI+Q energies to bicubic polynomials in 
tup and thh- Note that there is a substantial difference in both the saddle point 
geometry and in the zero-point correction when a more accurate bicubic polynomial 
is used instead of a simple biquadratic form to represent the surface. For example, 
at the CCI level, when the bicubic polynomial is used instead of a biquadratic 
one, the saddle point shifts to shorter tup by about 0.03 do and the zero-point 
correction increases by 0.09 kcal/mole. These differences are slightly larger for the 
CCI+Q case. These relatively large changes are a consequence of the extremely flat 
potential surface. Using the bicubic polynomial fits, we then determined the MEP 
as a function of thf by interpolation, and computed the normal mode vibrational 
frequencies at points along the MEP by standard methods [30]. The adiabatic 
ground-state potential curves for F+H 2 and F+D 2 are shown in Fig. 3. These 
curves were determined by adding the differential vibrational zero-point energy 
correction to the CCI(CCI+Q) energy along the MEP. The zero-point correction 
was taken as half the H 2 F(D 2 F) symmetric stretch mode frequency plus twice the 
H 2 F(D 2 F) zeroth-bending level, less half w e for The magnitude of the 

anharmonic correction to the symmetric stretch normal mode was estimated by 
fitting a Morse potential to the second and third derivatives of the potential energy 
profile along the direction of the symmetric stretch normal mode. In all cases, 
the value of u e x e obtained was within 20 cm -1 of the value for H 2 or D 2 , and 
thus including anharmonic effects for reactants and transition state would alter 
the zero-point energy correction by 0.02 kcal/mole or less. Since the shift in Tfjp 
between the classical and adiabatic barrier is approximately an order of magnitude 
greater than the corresponding shift in r # h , we have used internal coordinates in the 
present study. The differences between using internal coordinates and mass scaled 
cartesian coordinates for determining the reaction coordinate and symmetric stretch 
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directions (and their appropriate reduced masses) are expected to be insignificant 
owing to the small reaction path curvature. 

In Table IV we summarize the computed zero-point and tunneling contributions 
to the computed threshold for the F + H 2 and F + D 2 reactions. The results la- 
beled “E barrier” are for the analysis of the classical barrier location, and manifest 
a substantial zero-point effect. For CCI + Q, the saddle point symmetric stretch 
frequency (predominantly H 2 stretch) is 3768 cm -1 compared to 4401 cm -1 for 
free H 2 , while the lowest eigenvalue for the bend potential is 45.9 cm -1 leading to 
a zero-point correction of 0.64 kcal/mole. A one-dimensional tunneling correction 
was determined by fitting the classical potential Vj^jep(s) to an Eckart barrier [31] 
(we matched the curvature at the barrier of the MEP, and adjusted the exother- 
micity to minimize the root-mean-square deviation between the Eckart potential 
and Vmep( s ))- We assumed the experimental threshold energy corresponds to 
a reaction cross section equal to 25% of its value when the translation^ energy 
equals the barrier height. This tunneling correction further reduces the barrier by 
about 0.5 kcal/mole leading to a translation energy threshold of 1.0 kcal/mole for 
a 2.1 kcal/mole high barrier. As the experimental threshold is w 0.7 kcal/mole, 
this would suggest that the true barrier height is about 1.8 kcal/mole. However, 
when the same analysis is carried out on the adiabatic ground-state potential curve 
(labeled “adiabatic barrier”), the dynamical bottleneck is found to occur at r hf 
s= 3.16 and thh =1-42, a shift in rjjp of 0.25 00 from the computed saddle point 
geometry. At this geometry, the H 2 stretching frequency is 4178 cm -1 and the 
zero-point correction, while still lowering the barrier, is only 0.06 kcal/mole. The 
resulting adiabatic barrier height is 1.86 kcal/mole. Inclusion of an equivalent one- 
dimensional tunneling correction leads to an estimated translational energy thresh- 
old of 1.5 kcal/mole, which is about 0.6 kcal/mole lower than the classical barrier 
height. This implies that the classical barrier on the exact ab initio potential energy 
surface should be approximately 1.3 kcal/mole. This barrier height is significantly 
below the 2.35 kcal/mole estimate of Schaefer [3,4], and also slightly below our best 
estimate based upon direct calculation. 

We applied the same analysis to both the F + H 2 and F + D 2 reactions using the 
CCI and CCI+Q results in both cases. Interestingly, for F + D 2 the tunneling factor 
is in general smaller, but the shift in energy between the adiabatic barrier and saddle 
point is larger. The resulting translational energy threshold predictions, which are 
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nearly the same for both reactions, are estimated to be 0.67± 0.05 kcal/mole below 
the saddle point energy on the potential energy surface. On the other hand, if the 
zero-point and tunneling corrections are applied at the classical barrier geometry, 
the F+D 2 reaction has a 0.24 kcal/mole larger threshold energy. 

It is difficult to assess the uncertainty in this method of predicting the bar- 
rier height. Neumark tt al. [5] easily observed reactive scattering in the F + 
H 2 /D 2 crossed molecular beam experiments at reactive collision energies of 0.7- 
0.8 kcal/mole, but did not give an estimate of the magnitude of the cross section. 
We have chosen to assign the tunneling correction to the threshold as the trans- 
lational energy for which the tunneling probability is 0.125 (25% of the value at 
the barrier energy). While Steckler et al. [6] concluded that the Eckart tunnel- 
ing correction worked well for this system, there probably still is an uncertainty of 
±0.2 kcal/mole in this value (due in part to determining the magnitude of the cross 
section at the experimental threshold). The tunneling correction is significantly 
smaller for the adiabatic barrier than for the classical barrier, owing to a broader 
barrier with reduced curvature (as manifested by a smaller imaginary asymmetric 
stretch frequency). This effect is more pronounced for the CCI+Q results, which 
have significantly lower barriers than do the CCI. In addition, the change in cur- 
vature at the adiabatic barrier between F-I-H 2 and F+D 2 nearly counteracts the 
mass effect and results in approximately equal tunneling corrections for the two 
reactions. 

The other uncertainties in the analysis are easier to assess. The vibrational 
zero-point correction is sensitive to the saddle-point location, which should be ac- 
curately determined at this level of theory. One measure of the remaining error 
is the approximately 10% difference between the MRCI+Q and the CCI+Q force 
constants. For the vibrational frequencies along the minimum energy path we have 
used the energy levels of the actual CCI bending potential and demonstrated the 
unimportance of correcting for anharmonicity in the stretching coordinate. The 
energy difference between the adiabatic barrier and the saddle point is estimated to 
contain an uncertainty of ±0.1 kcal/mole. Combining these two error estimates, we 
believe the energy barrier in the collinear ab initio potential energy surface should 
be 1.0— 1.5 kcal/mole. 
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V. Conclusions 

Calculations in support of a classical barrier height of 1.4±0.4 kcal/mole are 
presented. The bending potential near the saddle- point geometry is found to be 
very flat and to possess a slight minimum at noncollinear geometries. Since our best 
bending potential is in relatively good agreement with the SEC results of Schwenke 
et a l. [8], it would be worthwhile to carry out quantum scattering calculations 
with this characteristic in the potential energy surface. A controversial point in 
the present theoretical calculations is the validity of applying a +Q correction for 
higher excitations to the CASSCF /MRCI calculations. Although this correction 
has no rigorous justification and is clearly an overcorrection when all important 
correlation effects can be included in the MRCI reference space, we believe that the 
+Q results are superior in this case and may well be an underestimate of the effects 
of higher excitations in the 9-electron treatments. This is based on FCI calibrations 
and in part on the belief that the full effect of the 2s correlation contribution to the 
barrier is only obtained at a very high level of correlation treatment. 

The contention that the barrier is lower than previous theoretical calculations is 
also strongly supported by variational transition state theory calculations employing 
a surface computed with the CCI method. In agreement with previous work [6], 
the true bottleneck for the reaction is shifted considerably out into the entrance 
channel region as compared with the saddle-point geometry, where the combined 
zero-point and tunneling corrections are relatively small. Using the experimental 
thresholds [5] to reposition the computed barrier height favors a result in the lower 
half of the range 1.4±0.4 kcal/mole. Therefore, our work argues in favor of previous 
studies [7-9] that have recommended a relatively low barrier, and disagrees with the 
recent contention [3,4] that the barrier could not be less than 2.35 kcal/mole. 
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Table I. FCI calibration of the classical barrier height 0 . 


A. At the FCI saddle point 



barrier exothermicity 



FCI * * * 6 

4.50 

28.84 



MRCI(300) 

5.18 

28.57 



MRCI(300)+Q 

4.43 

29.12 



MRCI(322)(0.05) 

5.00 

29.12 



MRCI(322)(0.05)+Q 

4.32 

29.21 



MRCI(322) (0.025) 

4.73 

29.17 



MRCI(322) (0.025) +Q 

4.51 

28.80 



MRCI(322) (0.01) 

4.71 

29.19 



MRCI(322)(0.01)+Q 

4.54 

28.84 



MRCI(522) (0.025) 

4.55 

29.41 



MRCI (522) (0.025) +Q 

4.32 

29.31 



B. At optimized saddle-point geometry 0 





r(F-H) 

r(H-H) 

barrier 

exothermicity 

FCI 6 

2.761 

1.467 

4.50 

28.84 

CPF 6 

2.801 

1.467 

4.40 

26.47 

MRCI(300) 

2.740 

1.476 

5.16 

28.57 

MRCI(300)+Q 

2.795 

1.467 

4.42 

29.12 

MRCI(311) 6 

2.722 

1.473 

5.15 

28.68 

MRCI(311)+Q 6 

2.773 

1.464 

4.39 

29.21 

MRCI(322) (0.025) 

2.761 

1.474 

4.70 

29.17 

MRCI(322) (0.025) +Q 

2.755 

1.475 

4.49 

28.80 


0 Energies in kcal/mole and bond lengths in oo. All calculations are done using 

the [4s3pld/2slp] basis set and correlating 7 electrons. The barrier is referenced to 

F. . .H2(50ao), and the exothermicity is computed using HF. . .H(50ao). 

6 Taken from Ref. 25 (FCI). 

c Geometry optimizations were done using a biquadratic fit to a grid of nine points. 


19 



Table II. Spectroscopic constants for the H^E^) and HF( 1 E + ) systems, and basis 
set superposition errors. 


H 2 



r e(«o) 

De(eV) 

SDCI 

1.402 

4.747 

Expt° 

1.401 

4.747 


HF 



r e (a 0 ) 

D e (eV) 

MRCI(200) 

1.731 

5.975 

MRCI(200)+Q 

1.735 

6.057 

MRCI(222) 

1.735 

6.064 

MRCI(222)+Q 

1.736 

6.070 

Expt° 

1.733 

6.123 


Superposition errors in kcal/mole 6 

H 2 (F)-SDCI c 

0.01 


F(H 2 )-SDCI 

0.15(0.16) 


F(H 2 j— MRCI(122) 

0.15(0.16) 


f-(h 2 )-sdci 

0.58(0.65) 


F~ (H 2 )— MRCI(222) 

0.73(0.78) 


F(H 2 )— SDCI(2p) 

0.06(0.07) 


F~(H 2 )-SDCI(2p) 

0.42(0.44) 


F(H)-SDCI(STO) <i 

0.35 



“Ref. 27. 

fc The MRCI(300)+Q saddle point geometry is used, r(F-H)=2.921 and r(H- 
H)= 1.450. The values in parentheses include a -fQ correction. 

“The nomenclature X(Y) denotes that X is computed with the Y ghost basis. 
d The Slater basis set is taken from Ref. 2. 
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Table III. Theoretical studies of the classical saddle-point geometry and barrier for 
the F 4 -H 2 reaction. 


Basis 11 

Level of treatment 

saddle-point 

barrier * * 6 

exothermicity 1 

A 

MRCI(300) 

2.793 

1.465 

3.47 

30.35 

A 

MRCI(300)+Q 

2.921 

1.450 

2.29 

31.00 

Slater 

MRCI(300) 

2.80 

1.46 

3.24 


A 

MRCI(322) 

2.914 

1.451 

2.63 

31.61 

A 

MRCI(322)+Q 

2.950 

1.450 

1.66 

30.47 

A 

MRCI(322)(2p) c 

2.899 

1.455 

2.99 

33.96 

A 

MRCI(322)(2p)+Q c 

2.910 

1.456 

2.42 

33.42 

A 

SDCI+Q 

2.638 

1.470 

3.87 

29.40 

A 

CPF 

2.939 

1.451 

2.44 

28.85 

A 

CCI(322) 

d 

d 

2.79 

31.8 

A 

CCI(322)+Q 

d 

d 

2.02 

30.7 

A+H(/)' 

CCI(322) 

d 

d 

2.73 


A+H(/) 

CCI(322) +Q 

d 

d 

1.95 


A-F(s) 

CCI(322) 

2.879 

1.447 

2.89 


A-F(j) 

CCI(322)+Q 

Expt. 

2.909 

1.445 

2.14 

31.73 


tt This letter “A” denotes the [5s5p3d2 flgj4s3p2d] basis described in the text. The 

Slater basis is taken from Ref. 2. 

6 The barrier is referenced to F. . .H2(50oo) 5 and the exothermicity is computed 
using HF. . .H(50ao). 

c These are 7-electron treatments (i.e. 2s correlation is excluded). 
d The MRCI(300)-I-Q saddle point geometry is used, r(F-H)=2.921 and r(H- 
H)= 1.450. 

e Denotes that a function of this angular momentum type has been added. 
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Table IV. Zero-point and tunneling effects on the barrier height of the F -I- H 2 and 
F + D 2 reactions. 


F + H 2 surface 


Classical Barrier 

Adiabatic Barrier 


CCI 

CCI + Q 

CCI 

CCI + Q 

T HF, &0 

2.879 

2.909 

3.070 

3.155 

thh, &o 

1.447 

1.445 

1.425 

1.421 

Barrier, kcal/mole 

2.888 

2.143 

2.639 

1.860 

Sym. stretch, cm -1 

3706 

3768 

4074 

4178 

Bend, cm -1 

68.5 

45.9 

68.5 

45.9 

Asym. stretch 0 , cm -1 

692t 

605t 

530 * 

371t 

Zero-point correction 6 , kcal/mole 

-0.602 

-0.643 

-0.076 

-0.057 

E barrier -1- zero point, kcal/mole 

2.286 

1.500 

2.563 

1.803 

Tunneling correction, kcal/mole 

-0.54 

-0.47 

-0.42 

-0.29 

Threshold, kcal/mole 

1.75 

1.03 

2.14 

1.51 

F + D 2 surface 


Classical Barrier 

Adiabatic Barrier 


CCI 

CCI + Q 

CCI 

CCI + Q 

iHF) 

2.879 

2.909 

3.010 

3.075 

*HH > OO 

1.447 

1.445 

1.430 

1.427 

Barrier, kcal/mole 

2.888 

2.143 

2.761 

1.997 

Sym. stretch, cm -1 

2623 

2667 

2811 

2876 

Bend, cm -1 

37.7 

19.1 

37.7 

19.1 

Asym. stretch®, cm -1 

512t 

448i 

428t 

334* 

Zero-point correction 6 , kcal/mole 

-0.488 

-0.532 

-0.220 

-0.233 

E barrier + aero point, kcal/mole 

2.400 

1.611 

2.541 

1.764 

Tunneling correction, kcal/mole 

-0.40 

-0.35 

-0.34 

-0.26 

Threshold, kcal/mole 

2.00 

1.26 

2.20 

1.50 


“From the normal mode analysis at the classical barrier, and computed from the 
curvature along the Eckart potential at the adiabatic barrier. 

6 For H 2 (D 2 ) we used w e =4401(3116) cm -1 , respectively, from Ref. 26. 
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Figure Captions 


Figure 1: CCI bending potentials for F+H 2 , where 0=0 corresponds to collinear 
approach. The bending curves are denoted by the CASSCF active space in C, 
symmetry notation. The energy zero is arbitrary, but the curves are postioned 
relative to each other on the basis of total energy. 

Figure 2: The CCI(82)+Q bending potential with respect to the reactants asymp- 
tote showing the first few vibrational levels for F+H 2 (solid lines) and F+D 2 (dashed 
lines). Vibrational energies with respect to the bending potential at 0=0 are given 
in cm -1 . 

Figure 3: Eckart potentials fit to the CCI+Q data for the classical potential along 
the minimum energy path (MEP) and for the adiabatic ground-state potential en- 
ergy curves for F-f-H 2 and F+D 2 . Values of the reaction coordinate s and tjjf 
are both given on the abscissa. The symbols represent the corresponding curves 
computed using the bicubic fit to the ab initio potential energy grid. All curves are 
relative to the energy of the reactants. 
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